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ABSTRACT 


This  dissertation  presents  the  development  of  a 
finite  element  method  of  analysis  of  interaction  problems 
between  an  elastic  structure  and  a  moving  fluid.  The 
fundamental  principles  of  the  method  have  been  formulated 
for  the  case  of  finite  non-isothermal  deformations  of  the 
structure  and  unsteady  flow  of  a  compressible,  viscous, 
heat-conducting  fluid.  Details  of  the  method  and  numerical 
examples  have  been  worked  out  for  the  simpler,  and  yet 
practically  important  conditions  of  isothermal  deforma¬ 
tions  of  the  structure  and  compressible,  inviscid, 
isentropic  flow  of  the  fluid. 

The  finite  element  formulation  of  the  method 
represents  a  numerical  solution  of  the  variational 
equations  of  the  problem.  The  necessary  variational 
principles  have  been  derived  for  the  fluid,  for  the  struc¬ 
ture,  and  for  the  joint  fluid-structure  system.  Since 
the  governing  equations  of  the  fluid  do  not  form  a  poten¬ 
tial  operator,  the  corresponding  variational  functional 
has  been  constructed  in  terms  of  the  actual  and  adjoint 
variables.  Similar  variational  functionals  have  been 
obtained  the  structure  and  for  the  iteraction  problem. 
Consequently,  the  resulting  finite  element  algorithm 
belongs  to  the  class  of  "weighted  residuals”  methods. 

A  quadrilateral  isoparametric  element  has  been 
developed  and  used  for  the  fluid.  In  order  to  remove 


the  instability  effects  inherent  to  the  problems  with 
convective  terms,  the  "upwind  weighting"  procedure  has 
been  employed.  The  structure  can,  conceivably,  be 
represented  by  any  type  of  elements;  in  the  present 
examples,  a  compatible  membrane  element  has  been  used. 

The  numerical  examples  demonstrate,  first,  the 
application  of  the  method  for  the  analysis  of  transient 
problems  of  gas  dynamics  in  one  and  two  dimensions. 
Subsequently,  transient  motion  of  a  membrane  interacting 
with  a  two-dimensional  gas  flow,  and  the  resulting 
perturbations  in  the  gas  itself,  have  been  analyzed. 
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1.  INTRODUCTION 

The  goal  of  the  present  method  is  to  analyze  fluid 
structure  interactions  of  a  general  nature.  The  structure 
is  an  elastic  solid  of  arbitrary  configuration.  The  fluid 
is  modeled  as  a  compressible  linear  Newtonian  fluid.  The 
method  is  not  restricted  to  small  displacements  in  the 
structure,  nor  small  perturbations  in  the  fluid. 

The  equations  of  motion,  conservation  of  mass,  and 
conservation  of  energy  form  the  basis  for  the  method.  The 
variational  principle  is  stated  in  general  coordinates  but 
most  work  is  in  Cartesian  coordinates.  This  change  is  for 
convenience  only  and  does  not  represent  a  limitation  of 
the  method.  The  structural  equations  allow  large  displace¬ 
ments  and  thermal  strains.  The  fluid  is  formulated  in 
spatial  coordinates,  and  the  resulting  convected  gradients 
are  included.  The  fluid  is  not  restricted  to  being  at 
rest  initially. 

The  problem  solution  is  based  upon  the  finite 
element  method  for  both  fluid  and  structure.  Contributions 
have  been  made  to  the  use  of  upwind  weighting  to  stabi¬ 
lize  the  numerical  solution  of  the  fluid  dynamics  equa¬ 
tions.  The  ability  of  finite  elements  to  approximate 
complex  shapes  allows  the  analysis  of  very  general  con¬ 
figurations.  Further,  the  extremely  simple  handling 
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of  the  fluid  structure  interface  is  a  major  reason  for 
the  success  of  the  numerical  application  of  the  method. 

To  demonstrate  the  viability  of  the  method,  several 
numerical  examples  are  presented.  The  Newtonian  fluid 
model  is  reduced  to  an  ideal  gas  to  simplify  the 
numerical  procedure.  As  a  result,  the  interface  boundary 
condition  is  reduced  to  full  slip  or  continuity  of  the 
surface  normal  velocity.  The  examples  include  one¬ 
dimensional  gas  dynamics  and  vibration  of  an  elastic 
membrane. 

The  finite  element  method  is  instrumental  in  the 
effective  solution  of  the  arbitrary  problems  at  hand. 

The  present  method  will  allow  the  use  of  virtually  any 
structural  element  in  current  program  libraries,  thereby 
allowing  the  extension  of  previous  problems  to  include 
interaction  effects.  Although  the  numerical  examples  of 
this  work  are  problems  of  unsteady  motion,  the  functional 
and  the  discretization  process  are  directly  applicable  to 
other  problem  types.  Some  caution  is  indicated,  however, 
as  the  equations  of  the  problem  are  nonlinear  differential 
equations.  Despite  the  limited  success  of  the  present 
work,  much  remains  to  be  done  in  the  numerical  solution 
of  problems  of  this  nature. 

By  far  the  most  common  formulation  of  structural 
finite  elements  is  based  upon  the  displacements  and 
velocities.  In  contrast,  fluid  dynamics  problems  are  more 


frequently  based  upon  the  velocity  potential  or  the  stream 
function  and  vorticity.  As  examples  of  fluid  structure 
interaction  become  more  common,  a  shift  to  the  use  of 
velocity  has  occurred.  The  use  of  a  common  basis  for  both 
structure  and  fluid  as  is  done  here  can  simplify  the 
mechanization  of  many  aspects  of  the  problem  solution. 

The  handling  of  the  interface  in  the  present  work 
is  simplified  greatly  by  treating  this  region  much  as  a 
structural  method  handles  an  inhomogeneous  medium.  That 
is,  as  long  as  the  velocity  field  is  continuous  across 
the  surface,  no  other  conditions  need  be  explicitly 
represented.  When  using  finite  elements  in  practice, 
this  means  using  compatible  elements  on  the  interface. 
Further,  it  requires  coincidence  of  the  structural  and 
fluid  nodes.  As  a  result,  the  common  difficulties 
arising  from  differing  mesh  densities  do  not  occur. 

A  result  of  the  present  work  is  the  demonstration 
by  numerical  examples  of  a  new  unified  method  for  fluid 
structure  interaction  problems.  The  method  is  unified 
in  the  sense  that  the  entire  problem  domain  is  formulated 
and  solved  by  the  same  procedure.  Differences  in  treat¬ 
ment  of  the  structure  and  fluid  domains  arise  as  a 
result  of  numerical  stability  requirements  rather  than 
the  variational  formulation. 


2 .  BACKGROUND 


There  appears  to  be  no  published  examples  of  fluid 
structure  interaction  involving  a  compressible  viscous 
fluid  in  unsteady  motion.  It  seems  little  attention  has 
been  paid  to  formulating  the  solution  procedure  for  such 
a  fluid  until  recently  (Nakamura,  1977  and  Chung,  1978). 
By  far  the  most  developed  fluid  solution  procedures  are 
for  potential  flow,  acoustic  fluids,  and  creeping 
viscous  flow.  Whatever  the  reason  for  this,  it  is 
probably  not  unrelated  to  the  fact  that  these  three 
classes  of  fluid  problems  have  simple  variational  formula 
tions.  Due  to  the  well  established  position  of  these 
procedures,  they  are  the  dominant  formulations  employed 
in  interaction  problems. 

For  example,  some  of  the  specific  fluid  structure 
combinations  that  have  been  tried  are:  elastic  shells 
and  acoustic  fluid,  both  by  finite  differences 
(DiMaggio,  Bleich,  McCormick,  1978);  elasto-plastic 
shells  and  acoustic  fluid,  again  by  finite  differences 
(Nikolakopoulou,  DiMaggio,  1978);  elastic  beams  with 
acoustic  fluid  by  finite  elements  (Zarda,  1976) ; 
elastic  bodies  and  incompressible  viscous  fluids  by 
finite  elements  (Zarda,  Chien,  Skalak,  1977);  elastic 
structural  response  to  internal  explosion  with  an 
acoustic  fluid  (Jacobson,  Yamane,  Brues,  1977) . 

In  addition,  a  lot  of  work  has  been  done  to 


uncouple  the  problem  by  employing  a  suitable  approximation 
for  the  fluid  response.  Two  of  the  most  common  are  the 
doubly  asymptotic  approximation  and  supersonic  thin 
airfoil  theory.  Examples  of  the  former  are:  elastic 
structural  response  to  shock  (Ranlet,  DiMaggio,  Bleich, 
Baron,  1977);  elasto-plastic  structural  response  to  shock 
(Atkatsh,  Baron,  Bieniek,  1978);  elastic  structural  response 
(Everstine,  1976  and  Everstine,  Schroeder,  Marcus,  1975); 
transient  structural  response  (Geers,  1975). 

The  use  of  supersonic  thin  airfoil  theory  allows 
the  local  pressure  to  be  related  to  the  structural  slope 
and  velocity.  Examples  are:  elastic  panel  stability 
(Kornecki,  Dowell,  O'Brien,  1976);  panel  flutter  (Olson, 
1970  and  Dowell,  1970);  large  deflection  plate  flutter 
(Dowell,  1966);  Lyapunov  stability  analysis  of  panels 
(Wang,  1966) ;  stability  of  cylindrical  shell  limit  cycle 
oscillations  (Evensen,  Olson,  1968) . 

Several  other  fluid  solution  procedures  are 
commonly  used  in  iterative  solution  processes.  By 
separating  the  problem  at  each  iteration,  the  method  used 
to  calculate  the  fluid  response  is  almost  entirely 
divorced  from  the  structural  solution  method.  For  this 
reason  fluid  solution  techniques  will  be  considered 
independently  from  the  structural  method  in  the  following. 

Perhaps  the  fastest  growing  inviscid  procedure  is 
the  doublet  lattice  method,  although  other  methods  of 
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singularity  distribution  are  also  used.  Some  are: 
doublet  lattice  method  for  interfering  wings  and  bodies 
(Rodden,  Giesing,  Kalman,  1970  and  Giesing,  Kalman, 
Rodden,  1972);  vortex  lattice  for  lifting  surfaces 
(Pittman,  Dillon,  1977);  doublet  lattice  (Jordan,  1978); 
distributed  sources  and  doublets  in  a  paneling  method 
(Dusto,  Epton,  Johnson,  1978);  distributed  sources  for 
bluff  body  wakes  (Parkinson,  Jandali,  1970);  distributed 
doublets  with  cubic  spline  interpolation  functions 
(Gotta,  van  de  Vooren,  1972). 

One  other  method  is  finding  acceptance  for  certain 
linear  inviscid  problems.  Known  variously  as  the  Green's 
function  method  or  the  boundary  solution  procedure,  its 
chief  advantage  is  to  reduce  volume  problems  to  surface 
problems.  The  former  has  been  applied  to  potential  flow 
where  the  surface  is  modeled  by  flat  quadrilaterals 
(Morino,  1973)  .  The  latter  has  been  used  in  surface 
wave  problems  (Zienkiewicz ,  1977)  and  other  problems 
(Zienkiewicz ,  Kelly,  Bettes,  1977). 

For  viscous  fluids,  the  majority  of  the  work  has 
been  done  in  very  low  Reynolds  number  flows.  Here  a 
simple  variational  principle  leads  to  minimizing  the 
dissipation.  Some  examples  are:  Newtonian  and  special 
non-Newtonian  fluids  (Delleur,  Sooky,  1961);  Galerkin's 
method  with  the  use  of  Lagrangian  multiplier  for 
incompressibility  (Argyris,  Mareczek,  1974); 
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minimization  of  dissipation  with  finite  elements 
(Oden,  Somogoyi,  1969).  Alternate  formulations  for 
viscous  flows  include  the  stream  function  and  vorticity 
(Baker,  1974  and  Tuann,  Olson,  1976) .  Series  solutions 
have  been  applied,  for  example,  to  elastic  spheres  in 
close  fitting  tubes  (Tozeren,  Skalak,  1978) . 

When  studying  higher  Reynolds  number  flow,  the 
use  of  upwind  differencing  (Steele,  Barrett,  1978)  or 
upwind  weighting  (Christie,  Griffiths,  Mitchel,  Zienkiewicz, 
1976)  brings  numerical  stability.  This  technique  is 

being  mentioned  more  frequently  in  new  texts  (Chung,  1978) 

♦ 

but  few  numerical  examples  beyond  two-dimensional  tempera¬ 
ture  convection  (Zienkiewicz,  1977)  are  given. 


t- 
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3.  VARIATIONAL  FORMULATION  OF  THE  PROBLEM 
3.1  Existing  Variational  Principles 

Examination  of  the  variational  methods  of  fluid 
dynamics  reveals  that  there  is  no  established  common  varia¬ 
tional  principle  valid  for  compressible  viscous  fluids  in 
unsteady  motion.  In  fact,  Finlayson  (1972)  has  shown  that 
the  Navier  Stokes  equations  do  not  comprise  a  potential 
operator.  Thus,  there  is  no  variational  principle  in 
the  conventional  sense. 

In  the  context  of  the  present  search  for  a  general 
procedure,  the  existing  variational  principles  suffer 
from  one  of  two  difficulties.  Most  frequently,  the 
problem  must  be  restricted  or  simplified  to  obtain  a 
functional.  An  example  is  the  assumption  of  small 
velocities  and  incompressibility  in  creeping  viscous 
flow.  Secondly,  side  conditions  may  be  imposed  on  the 
trial  functions.  As  the  problems  become  more  complex, 
it  is  difficult  to  include  all  the  equations  in  the  func¬ 
tional  without  introducing  unknown  multipliers.  Such 
conditions  are  then  left  as  constraints  on  the  trial 
functions  and  can  complicate  the  numerical  solution 
procedure . 


For  inviscid  compressible  flow  of  a  perfect  gas, 
Serrin  (1959)  presents  a  functional  based  primarily  on 
kinetic  and  internal  energies.  The  continuity  equation 
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and  the  isentropic  condition  are  added  by  using  multi¬ 
pliers.  This  leads  to  a  formulation  based  upon  velocity, 
density,  entropy  and  two  multipliers.  A  later  modifica¬ 
tion,  attributed  to  C.C.  Lin  by  Finlayson,  includes  the 
material  coordinates  and  a  multiplier  to  insure  that  these 
coordinates  do  not  change  in  rotational  flow. 

The  most  common  inviscid  principles  are  based  upon 
the  velocity  potential.  For  steady  incompressible  flow 
the  equation  to  be  solved  is  Laplace’s  equation  in  the 
velocity  potential.  The  well-known  functional  involves 
the  sum  of  the  squares  of  the  spatial  gradients  of  the 
potential . 

A  functional  for  steady  compressible  isentropic 
flow  based  upon  the  velocity  potential  is  easily  con¬ 
structed  (Carey,  1975  and  Norrie ,  de  Vries ,  1973).  The 
functional  is 


t  t*u,.r 

I  (<*»)  =  /v  j  fQ  p  (t)  dtdV 


(3.1-1) 


where 


P(t)  =  [1  +  1~-  win  -  t)]Y  ' 


(3.1-2) 


The  Euler  equation  for  this  functional  is  the  steady 
continuity  equation  in  terms  of  the  potential  <j> , 


(p$#i) fi  =  0 


(3.1-3) 
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Carey  (1975)  suggests  the  solution  should  be  based  upon 
discretization  by  finite  elements  followed  by  a  perturba¬ 
tion  expansion  to  handle  the  required  root.  The  first 
term  in  such  an  expansion  leads  to  the  incompressible  flow 
equation,  and  successive  terms  give  compressibility 
corrections . 

A  quasi-variational  principle  for  the  compressible 
flow  velocity  potential  equation  has  been  given  by  Norrie 
and  deVries  (1973).  For  two-dimensional  flow  the  equation 
is 


(<fc,^  —  a^)4>»  +2<Ji,  6,  (J>.  +  (<J> ,  ^  —  a^ )  <j> ,  =0 

*'x  xx  x  y  'xy  y  Y'yv 

(3.1-4) 


where  the  local  speed  of  sound, 


2  .  y  —  1  2.1-7  ^  *  2, 

a  =  a»o  +  1~T~  q°°  +  ~T~  U,x  4>'y) 


(3.1-5) 


The  functional  is  written  as 


I(v)  = ;  [rk  +  (v,!  +  v,.2)  +-^5  (v,vv,j)  dv 


V  2 


x  y 


-  / s  (Hv  +  Qv  +  |  av2]  dS 


(3.1-6) 


where 

G  =  4> 


'xy 'y' 


1  3  3  r 

H  =  — y  ( 4> ,  n  +  4,  n  )  +  — *■  ( ,  n  +4>,  n  ) 
3a2  x  x  y  y  **  x  y  y  x 


y 

(3.1-7) 


The  functional  of  Eq.  3.1-6  has  as  its  Euler  equation 
the  preceding  differential  equation  if  variations  are 
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permitted  in  v  only  and,  after  variations,  v  is  set 
equal  to  <f>.  In  practice,  the  functions,  G,  H  and  a  are 
evaluated  at  a  previous  iteration  and  treated  as  known 
on  the  next  iteration. 

Tor  slow  viscous  incompressible  flow  the  functional 
has  been  given  by  Zarda,  Chien  and  Skalak  (1977). 

The  incompressibility  condition  is  included  by  using 
the  pressure  as  a  multiplier.  The  functional  is 


I(VP)  =  V^aB^g-PVct1  dV+/S2VatadS  (3a-8) 


where 


d«e  *  I  <va,e  *  vb,«)' 


(3.1-9) 


velocity  is  specified  on  S1  and  tQ  is  a  prescribed 
traction  on  S2-  The  Euler  equation  is 


2udct8,e  "  P,a  ° 


(3.1-10) 


while  the  boundary  condition  on  S2  is 


2yd  DnD  -  pn  =  t 
aB  B  c  a  a 


(3.1-11) 


Steady  heat  conduction  with  constant  thermal  con¬ 
ductivity  has  a  variational  principle  very  similar  to 
incompressible  potential  flow.  The  functional  is 


KT)  =  /v[kT,aT,a  -  Tr]  dV+  /^Tq^dS  (3.1-12) 


where  q  is  a  prescribed  heat  flux  on  surface  S_,  r  is 

8  4 c 

the  rate  of  heat  generation  and  T  is  prescribed  on  Sj^. 


...  \«l  WUlWSMlMUk.  . 


The  Euler  equation  is 


kT.  +  r  =  0  (3.1-13) 

do 

and  the  boundary  condition  on  S2  is 

'  kT'ana  *  13.1-14) 

The  unsteady  heat  conduction  equation  is  non- 
selfadjoint  due  to  the  presence  of  the  first  time 
derivative.  As  a  result  the  operator  is  not  potential 
and  no  variational  principle  exists.  With  certain 
restrictions,  however,  the  equation  can  be  transformed 
into  a  symmetric  form.  For  linear  thermoelasticity, 
for  example,  the  Laplace  transform  in  time  can  be  used 
to  render  the  equations  symmetric  (Nickell ,  Sackman,  1968). 
The  method  has  the  advantage  of  having  the  initial  con¬ 
ditions  contained  explicitly  in  the  functional. 

The  discussion  of  variational  principles  would  not 
be  complete  without  a  brief  mention  of  the  Principle  of 
Minimum  Total  Potential  Energy.  The  principle  is 
applicable  to  isothermal  motion  of  an  elastic  solid. 

It  states  that  minimizing  the  sum  of  the  strain  energy 
and  the  potential  of  external  forces  is  equivalent  to 
satisfying  the  equations  of  motion.  Use  of  the  principle 
requires  the  specification  of  the  strain  energy  func¬ 
tion.  For  isotropic  materials  with  nonlinear  con¬ 
stitutive  equations,  a  suitable  form  that  will 
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approximate  the  behavior  of  the  material  can  usually 
be  found  by  using  functions  of  the  strain  invariants. 


3.2  General  Method  of  Construction  of  Variational 
Functionals 

Several  of  the  more  well  developed  areas  of 
structural  and  fluid  dynamics  have  in  common  the  exis¬ 
tence  of  functionals  and  variational  principles.  Linear 
elasticity,  heat  conduction  in  solids,  slow  viscous  flow 
and  potential  flow  are  just  a  few  of  the  areas.  In  most 
cases,  determination  of  the  existence  of  the  functional 
for  a  given  differential  operator  has  been  a  matter  of 
trial  and  error.  If  the  operator  represents  the  Euler 
equations  of  a  functional,  that  operator  is  the  gradient 
of  the  functional  in  a  variational  sense.  The  operator 
is  then  said  to  be  potential.  Finlayson  (1972)  presents 
a  method  based  upon  Frechet  differentials  to  easily 
determine  whether  an  operator  is  potential .  He  uses 
the  method  to  demonstrate  that  a  functional  does  not 
exist  for  steady  motion  of  an  incompressible  viscous 
fluid  unless 

v  x  (V  x  v)  =  0  or  v  •  Vv  =  0  (3.2-1) 


where  v  is  the  velocity.  The  following  is  a  summary  of 
Finlayson's  work  leading  to  the  use  of  adjoint  variables 
to  construct  a  functional  for  non-potential  operators. 
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(3.2-7) 


The  Euler  equations  are 


6u  :  N (u)  -  f  =  0 


~  i  *  *  * 

6u  :  Nu  -g  =  N  (u,u  )  -  g  =  0 


(3.2-8a) 
( 3 • 2-8b) 


Thus, if  the  problem  is  expanded  to  include  the  adjoint 
equation,  a  variational  principle  exists. 

3.3  Compressible  Newtonian  Fluid 
Constitutive  Relations 

The  stresses  at  a  point,  P^g,  consist  of  a 
pressure,  p,  and  viscous  stresses,  t  0,  such  that 


PaB  =  ■  p<5a6  +  tae 


(3.3-1) 


The  Newtonian  fluid  model  for  the  viscous  stresses  assumes 


t  „  =  Xv  6  _  +  2yd  D 
aB  y ,  Y  aB  aB 


(3.3-2) 


where  X  and  y  are  coefficients  of  viscosity,  v^  is  the 
velocity  and  the  rate  of  strain  tensor  is 


daB  2  *Va»B+  V6,a^ 


(3.3-3) 


The  equation  of  state  for  the  pressure  can  be  written  as 


p  =  pRT 


(3.3-4) 


For  a  perfect  gas,  the  entropy  function  per  unit  mass  is 


S  =  — —  log 

Y  -  1  PT 

and  the  internal  energy,  per  unit  mass, 


E  =  CyT 


(3.3-5) 


(3.3-6) 


1 
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If  0^  is  the  specific  heat  at  constant  pressure,  and 
Cv  is  the  specific  heat  at  constant  volume,  R  and  y  are 


defined  as 


R  *  C  -  C 
P  v 


y  =  C  /C 
P  v 


The  heat  conduction  law  is  assumed  as 


q  =  -  kT, 

a 


(3.3-7) 


where  q  is  the  heat  flux  vector,  k  is  the  thermal 

a 

conductivity  and  T  is  temperature. 

The  Governing  Differential  Equations 

The  equation  of  motion  for  the  fluid  in  spatial 


coordinates  is  assumed  in  the  form 
Dv 


Dt  ~  Pa8 , 6  pfa  "  0 


(3.3-8) 


where  f&  is  the  body  force  vector  per  unit  mass.  After 
using  the  Newtonian  fluid  model  and  the  perfect  gas  law, 
the  equation  of  motion  becomes 


is 


Dv 

0  D#  -  >ve,Ba-2udt.B,B  +  (RT0,'tl-»£a  *°  (3-3-’> 


The  equation  of  continuity  or  conservation  of  mass 


81  +  °Va  *  0 


(3.3-10) 
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The  energy  equation  is 

,DE  1  , 

P  (Dt  2  Dt  (PaBVa)  '  6 


q  +  pf  v 
^a,a  a  a 


(3.3-11) 


Dotting  the  equation  of  motion  into  the  velocity  and  sub¬ 
tracting  gives  the  thermodynamic  energy  equation. 


DE  _ 

0  Dt  =  P-«V-  *  ‘  <J. 


a8  a(  £  ^a,a 


(3.3-12) 


After  inserting  the  Newtonian  fluid  model,  the  perfect 
gas  law,  and  the  heat  conduction  law,  the  result  is  the 
temperature  equation. 


C  +  RTv  -  X 

v  Dt  a , a 


Va,aV8,B  *  2pdaBdaB 


k  T 


aa 


=  0 

(3.3-13) 


The  Functional 

In  the  following,  the  form  of  the  functional  in 
Eq.  3.2  -7  corresponding  to  the  preceding  equations 
is  given.  The  boundary  conditions  are  extracted  from 
the  volume  integral  by  integration  by  parts  of  the 
stresses  and  the  temperature.  In  addition,  adjoint 
boundary  functions  are  added  to  enhance  the  appearance 
of  symmetry.  Let  V_  denote  the  fluid  volume  and  S  its 

r 

surface.  S-,  is  that  part  of  S  where  tractions  t  are 
2  * 
specified  and  Sp  that  part  where  the  heat  flux 

is  specified. 


V. 


18 


The  functional  is  given  by 

*2  *  Dv  .  *  * 

I  =  /  /  {v  r  -~  +  Av  vD  e  +  2pd  0d  -  RTpv 

...  a  Dt  ot,a  B,B  aB  <*B  a. a 

Z1  VF 

*  *  Do  *  *  DT  * 

"  vaDfa  +  P  Dt  +  p  pvct,a  +  TCv  “+T  RTva,a 

“  T‘>,Va,aVe,e-T*2u'da6dae  +  k,T'«T'a}dVdt  " 


u  9  *  *  A  *t  ^ 

/  /  {v  t  +  v  t  }  dSdt  + 

.  c  a  a  a  a 

C1  °p 
i  1 2 


+  I  2  J  {T*q  n  +  Tq*n  }  dSdt 

.  a  a 

i  c  2 


(3.3-14) 


The  Euler  Equations  and  Boundary  Conditions 
For  the  functional  defined  in  Eq.  3.3-14,  the 
Euler  equations  in  the  primary  variables  are 


Dv 


’V  CD r"  XvB,6a-2l,aa6,6  +  (RT<»  m  0 

(3.3-15a) 


6  p  :  +  p  v  =0 

Dt  a,  a 


(3. 3-15b) 


6T* :  C  —  +  RTv  -A  V  v_  B  -  2y  d  Rd  R  -  k  T,  *  0 
v  Dt  a, a  a, a  B,B  aB  aB  cxa 

(3.3-lSc) 


The  remaining  Euler  equations,  mixed  in  the  primary 
variables  and  the  adjoint  variables,  are 
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K:  ♦  (v«'pv,),, ,  +  V  <* ,  ♦  fij  + 

*  f*T*T)  -  «7t-V,,>;<<  -  •*y(r’J„),t  - 


-  w« 


Vp  v«  'T,£v% 


( 3 . 3-15d) 


{p  :  P*  -  *  v/i  *  rt^  +  v,  p*  -  o 

JT:  cX  *  (T\V;  *  fi  P  C  '  *  rV-.«  + 


(3. 3-15e) 


+  kj  9  0 


( 3 . 3-15f ) 


The  boundary  conditions  are 


07,  SFi:  X  y, +  2^7),  -RTp«,  =  £  (3.3-lta) 

*  tew  +  p*p”«  f 
+  t7rT«.  -  =  *« 


»  J. : 


-  k’X  a, 

'*  *< 


J  77 

I*  * 


-  -  TXtV  =? 


-  ?*» 


( 3 . 3-16b) 
(3.3-16c) 

( 3 . 3-16d) 
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3.4  Inviscid  Isentropic  Fluid 

Assumptions  and  Constitutive  Equations 

This  fluid  model  is  a  special  case  of  the  compress¬ 
ible  Newtonian  fluid  discussed  in  the  preceding  section. 
Since  this  model  is  used  frequently  in  various  analyses, 
its  assumptions,  equations  and  functional  will  be  written 
here  explicitly. 

The  assumption  of  negligible  viscosity  is  expressed 
by  setting  the  constants  >  and  u  of  the  Newtonian  fluid 
model  equal  to  zero.  Consequently,  the  viscous  stress 
tensor  t  c  vanishes  identically. 

The  equation  of  state  for  a  perfect  gas  is 


p  =  CRT 


(3.4-1) 


With  the  assumption  of  perfect  behavior  the  entropy 
function  is  given  by 


S(P,P)  =  log 


(3.4-2) 


Again,  0^  is  the  specific  heat  at  constant  pressure,  Cv 
is  the  specific  heat  at  constant  volume  and 


R  =  C 


(3.4-3) 


y  =  C  /C 
P  v 


(3.4-4) 


The  isentropic  assumption  for  a  moving  fluid  takes  the  form 


DS 

Dt 


0 


(3.4-5) 
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Thus  the  energy  equation  in  terms  of  entropy 


T 


DS 

Dt 


1 

P 


q 


a ,  a 


(3.4-6) 


leads  to  q  =  0  or  no  heat  conduction. 

,  a 

Body  forces  will  not  be  considered,  but  may  easily 
be  included. 

Full  slip  will  be  assumed  at  all  boundaries,  and 
only  the  normal  velocity  of  the  fluid  will  be  prescribed. 

It  will  be  assumed  that  the  internal  energy  per 
unit  mass,  E,  is  a  function  of  temperature  only,  in  the 
form 

E  =  C  T  (3.4-7) 

v 


Governing  Differential  Equations 

Without  viscous  stresses,  the  equations  of  motion 
reduce  to 

Dv 

p  'DF  +  (RTp)  'a  =  0  (3.4-8) 

when  the  perfect  gas  law  is  used  to  eliminate  the 
pressure. 

The  continuity  equation  is 


^  +  pv  =0 
Dt  a ,  a 


(3.4-9) 


With  the  isentropic  assumption  and  the  lack  of  dis¬ 
sipation,  the  temperature  equation  is  simply 
DT 

C  kt  +  RTv 
v  Dt  a, a 


0 


(3.4-10) 
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The  Functional 


Insertion  of  the  preceding  equations  into  Eq.  3.2-7 
leads  to  the  functional  to  be  presented  in  the  following. 
Without  the  viscous  stresses  and  heat  conduction, 
integration  by  parts  is  performed  only  on  the  pressure 
gradient.  Thus,  the  only  surface  integral  in  the 
functional  contains  a  specified  pressure. 

The  inviscid  isentropic  flow  functional  is 
it 


-  rtpv-  +  p'2 


*  e' pe¬ 


er 


+  fcv» 7  *  r'RTv^J  JvJi  - 


-f  X  L  e' t**  *  e  J  Js  J* 


(3.4-11) 


The  Euler  Equations  and  Boundary  Conditions 
For  the  functional  defined  by  Eq.  3.4-11,  the 
Euler  equations  on  the  primary  variables  are: 


: 

Dv* 

P  Di 

+ 

O 

(3 . 4-12a) 

££. 

Di 

f 

pe.,  * 

o 

( 3 . 4-12b) 

cr: 

C  - 
^  Oi 

+ 

0 

(3.4-12C) 
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The  Euler  equations,  mixed  in  the  primary  and  adjoint 


variables 

,  are 

Sv„: 

PMV 

Z>* 

-pv>»,a  +  pv/v,,«  - 

+  T%  Td  -CfR  TJV  *o 

( 3 . 4-12d) 

SP: 

_V 

Di 

Dy„ 

*  v*1  1  -  ptV  so 

+  v*  ^ 

( 3 . 4-12e) 

ST: 

T)T* 

^Di 

-CA.,*  RTV,  -Rpv4'° 

( 3 . 4-12f ) 

The  boundary  conditions  are 

OH  ( 3 •  4-13a) 

£pVari e  +T'RTnd  4  p#p*nH  «•  p\  (3.4-i3b) 
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3.5  Elastic  Solid 

The  equations  for  thermodynamics  of  an  elastic 
solid  are  presented  here  only  for  completeness.  The 
functional  given  in  this  section  is  not  that  of  the 
Principle  of  Minimum  Total  Potential  nergy  but  instead 
is  the  adjoint  functional  constructed  by  the  present 
method.  An  example  will  be  given  in  a  subsequent 
chapter  of  an  elastic  membrane  with  geometric  non- 
linearities.  It  will  be  shown  that,  in  this  case,  the 
results  are  identical  to  those  of  the  Principle. 

The  Constitutive  Equations 

Let  U  be  the  internal  energy  per  unit  mass,  T  the 
temperature,  and  S,  the  entropy  per  unit  mass.  The  free 
energy  A  is  defined  as 

A  =  U  -  ST  (3.5-1) 


For  a  linear  elastic  material. 


°0A  *  7*  leaa’2  *  tJeabeab 


the  free  energy  will  be 

-  -  k  Yt2  (3.5-2) 

aa  2 


where  A  and  u  are  the  Lam6  constants,  B  is  the  pressure 
coefficient  and  y  is  a  specific  heat  parameter.  Thus, 
the  material  stress  tensor  is  found  to  be 


.  _  0  3  A 

sab  p0  3  e 


ab 


-  Xecc6ab  +  2l,eab  '  6T6ab 


(3.5-3) 


Further,  the  entropy  function  is  given  by 


3_A 
3  T 


+ 


1_ 

PO 


YT 


S 


(3.5-4) 
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In  the  above,  the  material  strain  tensor  is 


(u  ,  +  u,  +  u  u  ,  ) 
a,b  b,a  c,a  c,b 


(3.5-5) 


The  heat  conduction  law  is  assumed  in  the  form 


where  qQ  is  the  heat  flux  vector,  and  k  is  the  thermal 

cl 

conductivity. 

Governing  Differential  Equations 
The  Equation  of  Motion 

The  equation  of  motion  for  the  elastic  solid  in 
material  coordinates  is  assumed  in  the  form 

Uab(Scb  +  uc,b> ' 'a  +  °0fc  ‘  p0ac  *  0  (3.5-7) 

where  f  is  the  body  force  vector  per  unit  mass,  and  a 
c  c 

is  the  acceleration  vector. 

The  Temperature  Equation 

The  local  form  of  the  First  Law  of  Thermodynamics 
minus  the  rate  of  change  of  mechanical  energy,  in 
material  coordinates,  is 

00°  -  sabaab  +  ^  *  c0r  *  0  (3-5-8) 

a, a 


The  rate  of  change  of  internal  energy  U,  can  be 
expressed  in  terms  of  the  free  energy  A,  entropy  S  and 
temperature  T,  as  follows: 


,  Ik- 


•  •  JA 

U  *  A  +  ST  +  ST  =  jfT  +£-e#l  +  ST  *ST 


(3.5-9) 


With  this  substitution  for  u,  the  energy  equation  becomes 


p0st  +  q0 


a  ,a 


-  pQr  =  0 


(3.5-10) 


For  the  specific  free  energy  function  given,  the  entropy 
rate  is 


S  = 


(3.5-11) 


and  the  equation  takes  the  form 


yTT  -kT,  -  pnr  +  BTe  =0  (3.5-12) 

aa  u  aa 

The  Functional 

The  functional  presented  here  is  obtained  from 

Eq.  3.2-7  by  insertion  of  the  preceding  equations.  It 

differs  from  the  Principle  of  Minimum  Total  Potential 

Energy  primarily  in  the  inclusion  of  the  temperature 

equation.  Integration  by  parts  is  performed  on  the 

acceleration,  the  strain  rate  term  in  the  temperature 

equation,  the  stresses  and  the  heat  flux  gradient.  The 

latter  two  will  contribute  to  the  boundary  conditions 

on  Sc  and  S_  ,  respectively,  after  the  variation.  The 
S2  S2 

specified  adjoint  boundary  functions  are  included  to 
enhance  the  appearance  of  symmetry. 
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The  functional  I  is  defined  by 

J  +  V  +  ‘Se/c  +  e.u, c*\+ 

+  klCTi  *  rT*TT  -  p.T'r  -  4- 

*{\l<K  *  + 

+  //  f  T^.  4  TfcttJjSM  (3.5-13. 

*■>  ■«st_ 


Euler  Equations  and  Boundary  Conditions 
For  the  functional  defined  in  Eq.  3.5-13,  the 
Euler  equations  in  the  primary  variables  are 

K'LS*Al  *  ♦  fti  -  &“<:  (3.5-14a) 

Sf:  rTT  -kX^-fir  *  Qe.iT  (3.5-14M 


The  Euler  equations,  mixed  in  primary  and  adjoint 
variables,  are 

S<4  ■  %k  4  (ic(  +  fiSfci  (Hi*»  +  + 

V“i*  •*  < +  <c«m+  v  <t)  - 
-e(Tr,)'JJh  -  *t*  =  °  ,3-5 


ST :  -ytt*  -*T#  +  ee^T*  *  o  o.s-Md) 


The  boundary  equations  are 


OX  SSi:  S^a.k  +  Ue,h)\  =  (3.5-15a) 

Me'A‘r'‘ +(V  *\AxiuK+  + 
+/'K<- +  “A +  <tV  *  v )  - 


*  o 

ll 

(3.5-15b) 

>\ 

o 

n 

(3.5-lSc) 

b 

4. 

TV 

( 3 . 5-15d) 

h 

4. 

3.6  Variational  Principle  for  Fluid  Structure  Interaction 

The  functional  for  fluid  structure  interaction  to  be 
presented  in  this  section  is  constructed  from  the  governing 
differential  equations  of  each  region.  The  condition  of 
continuity  on  the  interface  is  expressed  as  an  admissi¬ 
bility  requirement  on  the  trial  functions  and,  therefore, 
does  not  appear  in  the  functional. 

In  the  following,  general  coordinates  have  been  used 
with  the  vertical  bar  signifying  covariant  differentiation, 
and  the  metric  tensor  denoted  by  g^  or  The  use 

of  general  coordinates  in  this  principle  is  a  simple 
extension  of  the  previous  sections  and  is  justified 
solely  on  the  grounds  that  it  contributes  greatly  to  the 
generality  of  the  formulation. 
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Consider  the  region  R  divided  into  two  volumes,  Vp  and 

Vg  and  bounded  by  surface  S  (Fig.  3.1).  Let  Vgbe  that  part 

of  the  region  containing  the  structure.  The  portion  of 

S  bounding  Vg  is  Sg.  The  remaining  part  of  the  region, 

Vp,  contains  the  fluid.  The  portion  of  S  bounding  Vp 

is  Sp.  The  surface  S  will  be  partitioned  two  ways: 

first,  S  ,SC  ,S„  ,S„  ;  and  second,  Sc  ,S_  S_,  ,S_  . 

S1  S2  F1  F2  S1  S2 '  F1  F2 
Each  partition  will  account  for  the  entire  surface  S. 

Make  stationary,  the  functional 


1  =  /  /  [-<i.  \  +  P.“'  + 

v5 

.  LT1^  f^vTT  _ _ T  *  _ 


+  +  r  ITT  -  P.rT  - 

-T*f  + 

F 

*  +  +  p'p + 

+  T*CV~  +  RT'Tv^  -  JT/,,  v',,  - 

*  *'T*T + 
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-  f  l  [  < r  +  ^P]JSU  + 

*,  Vi 

+  /‘X  fr'tf  t,*  +Tf*\]jsji 

6f  a c 


c  *  * 

among  functions  u  ,uc#T,T  ,  continuous  and  defined  on 
Vg  and  Sg,  satisfying 


C  A  0  ^  A  A 

u  =  u  ,  u„  =  u„  on  S.  ( 3 . 6-2a) 

c  c 

T  =  T  ,  T*=  T*  on  Sc  (3.6-2b) 

bl 


,  n  *  *  * 

and  among  functions  v  ,va,o,c  ,T,T  ,  continuous  and 
defined  on  Vp  and  Sp,  satisfying 


v 


Q 


T  = 


* 

T  =  T 


where 


on  S 

F1 

on  S 

'  1 


( 3 . 6  -  3a ) 
(3.6-3  b) 


D*  3  •  B 
Dt  3  t  v  ’  |  B 


(3.6-4) 


Further,  let  the  interface  Sj.  be  the  common  boundary 
between  Vg  and  Vp.  Admissible  structural  velocities  and 
fluid  velocities  must  be  continuous  across  Sj,  and 
admissible  fluid  temperatures  and  structural  temperatures 
must  be  continuous  across  S j . 
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Stationarity  of  this  functional  is  equivalent  to 
the  satisfaction  of  the  governing  differential  equations 
and  assures  continuity  of  the  traction  vector  across  S^. 
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4.  NUMERICAL  METHOD  FOR  INVISCID  COMPRESSIBLE  FLUID 
AND  AN  ELASTIC  SOLID 

4.1  Spatial  Discretization  of  the  Fluid  Functional 
Matrix  Equations 

Proceeding  by  the  finite  element  method,  the  fluid 

volume,  V_ ,  is  subdivided  into  elements  of  volume  V  , 

'  F  Fe 

such  that  there  is  no  gapping  or  overlapping,  and  the 

entire  volume  is  included.  The  functional  I  will  be 

taken  as  the  sum  over  all  elements  of  the  element 

functionals  I  . 

e 

Within  each  element,  the  variables  will  be  approxi¬ 
mated  by  the  product  of  interpolation  functions  and 
nodal  parameters.  These  parameters  will  be  taken  to 
be  the  value  of  the  variable  at  the  spatial  point  or 
node  of  the  element. 


Specifically,  the  approximations  are 

*  *  * 

v  =N.v  v  =N.v 

a  1  ai  a  1  at i 

*  *  * 

P  =  L .  p .  p  =  L .  p . 

11  11 

*  *  * 

T  =  M.T.  T  =  M.T. 

11  li 


(4.1-1) 


Notice  that  the  following  two  expressions  are  equivalent: 


*  *  * 

v  =  N  .  v 
a  1  ai 


*  ★  ★ 

V  =  [N  ]  {v  } 
a 


The  vector  of  nodal  velocity  parameters  {v)  will  be 

ordered  in  coordinate  sequence  by  node.  Let  [n  1 

xy 

be  a  row  matrix  such  that 


i 
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i 


v  =  [N  1  (vl  (4.1-2) 

a,  a  xy 

Then  the  elemental  functional  becomes 

i=//  + 


2,  V,1 

£  ^ 
cl  . 


+//  lr,fKJT[M][ti)ivji  + 


4 

4 


f  /4 A  fr*f  - 

A  v 


-//  nvfw 

*»  s£ 


(4.1-3) 


Due  to  the  linearity  of  the  functional  in  the  adjoint 
variables,  enforcing  stationarity  with  respect  to  the  adjoint 
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H  =/  {MfM1:  [LjMMUiMpy 

e  %  *  J  (4.1-6f> 

(4 . l-6g) 

The  Bilinear  Quadrilateral  Fluid  Element 

The  isoparametric  quadrilateral  with  linear  inter¬ 
polation  functions  was  chosen  since  it  is  the  lowest 
order  quadrilateral  that  will  produce  a  continuous 
variable  field.  The  element  matrices  are  evaluated  by 
numerical  integration  using  a  two  point  Gaussian  quadra¬ 
ture  in  each  coordinate.  The  mapping  functions  and  all 
variable  interpolaton  functions  were  chosen  to  be 

Ni*7(l  +  ^i)  (1  +  nrii)  i  -  1,2, 3,4  (4.1.7) 

The  end  product  of  the  element  routine  is  an  elemental 
contribution  to  the  global  rate  vector.  Given  the  current 
value  of  the  nodal  parameters,  the  element  contributions 
are  calculated  and  summed  by  connectivity.  Lumped  masses 
are  used  to  simplify  the  procedure. 

4.2  Upwind  Weighting 

Convected  Inertia  Terms 

The  recent  work  of  Heinrich,  Huyakorn,  Mitchell 
and  Zienkiewicz  (1977)  forms  the  basis  for  the  present 


use  of  upwind  weighting  in  finite  elements.  The  authors 
have  established  the  nature  of  the  numerical  instability 
associated  with  one-dimensional  transport  problems.  By 
employing  an  asymmetric  or  "upwind  weighted  element," 
they  were  able  to  demonstrate  the  values  of  the  winding 
parameter  that  lead  to  unstable  solutions  to  the  difference 
equations.  Their  work  indicates  that  large  coefficients 
of  the  gradients  are  the  source  of  the  instability. 

The  equation  solution  is  stabilized  by  modifying 
the  weighting  functions  of  a  weighted  residual  method. 

For  example,  for  a  one-dimensional  element,  the  weighting 
functions  are 

W.  =  W^(x,a)  =  N^(x)  +  aF(x)  (4.2-1) 

The  authors  chose 

F  (x)  =  -  -4  x(x-h)  (4.2-2) 

h 

where  h  is  the  mesh  size  (Fig.  4.1).  The  effect  of  the 
upwinding  is  to  shift  the  elemental  emphasis  to  the 
direction  of  upwind.  This  is  accomplished  by  giving  a 
the  sign  of  the  velocity  v  along  the  element.  Their 
work  shows  that  a  certain  value  of  a  will  lead  to  minimum 
error  in  the  difference  equations.  A  sizable  reduction 
in  error  is  reported  when  using  an  optimal  value  of 
upwinding  instead  of  full  upwinding. 

The  authors  have  carried  their  numerical  work  to 
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two  dimensions  although  it  was  not  possible  to  extend  the 
analytical  work.  Since  two-dimensional  interpolation 
functions  are  formed  from  products  of  one-dimensional 
interpolation  functions,  the  two-diifensional  weighting 
functions  are  formed  similarly.  Introducing  6  as  an  n 
direction  winding  parameter,  the  functions  for  the 
bilinear  isoparametric  quadrilateral  are 

N;  *  i  ( l  +  C Ci  )  i  ( 1  *  tty)  (4.2-3) 

wt  iHtMfrhfaO+fXi-z) (4.2-4) 

Notice  that  these  weighting  functions  reduce  to  the 
interpolation  functions  for  the  case  of  no  winding. 

The  authors  note  that,  in  a  true  method  of  weighted 
residuals  fashion,  these  winded  weighting  functions  are 
applied  to  all  terms  in  the  residual.  The  analogous 
development  foar  the  present  set  of  nonlinear  equations 
has  not  been  demonstrated. 

Typical  Nodal  Equation  in  One  Dimension 
To  see  the  effects  of  upwind  weighting  when  used 
in  the  present  set  of  nonlinear  equations,  the  typical 
nodal  equation  for  a  one-dimensional  assemblage  of 
elements  will  be  presented.  The  elements  are  line 
elements  of  unit  length  using  linear  interpolation  func¬ 
tions.  The  upwind  weighting  functions  will  be  applied 
to  the  convected  gradient  terms  only. 


rtfn  rwMTt.-a-.g  mcar~.  - rsmj rr. 
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Consider  two  line  elements  of  unit  length.  Let  the 
common  node  be  labeled  2,  the  upwind  node  be  1,  and  the 
downwind  node  be  3.  Each  node  will  have  three  nodal 
parameters:  u^,p^,T^.  By  straightforward  evaluation  of 
the  matrix  equations,  the  three  typical  nodal  equations 
are  found  to  be 


+  (£{>,+  *■  + 

+ K-v.te  p .  ( 1 4*0  ♦  ,i  ej'  <  HJ  i  + 

<{(*-*, *?«>  *th(>  44 + * 

- : R  f i  ■ i;  p,  *  i  ‘ T,  ft  +  f T  ft  -  f H  ft  -  f  S  p,  - f t,  P, ) 1 = < 0 

(4 . 2-5a) 

+  3  OW«)0"i*()V} 

-(tc.*ihM  *  +  ^p« 

(4 . 2-5b) 


V 


1 


it  Tt  +  1%  +l(Tz-T, )(,  *§<)\(  +  fl(\-TiXn$><)  4- 

+  jfo -iiVi-fok  +  {n-r^i-hy,  - 
~rMr,  +cAT‘ -jt»K  ♦  J(* t'iW,-o 

(4.2-5C) 

For  a  local  velocity  positive,  or  from  node  1  to 
node  3,  a  =  1.0.  The  gradients  from  the  left  or  upwind 
element  receive  a  greater  weight  than  corresponding 
gradients  from  the  downwind  element.  In  the  momentum 
equation,  the  upwind  gradient,  v2~vi'  ^as  weight 
increased  by  a,  that  part  including  receiving  the 
largest.  The  downwind  gradient,  v^-v^  receives 
diminished  weight  when  multiplied  by  v2,  and  even  further 
reduced  weight  when  multiplied  by  v^.  In  both  terms, 

P2  receives  more  weight  than  p^.  For  negative  velocities, 
a  -  -  1.0  and  all  effects  are  reversed. 

Typical  Nodal  Equation  in  Two  Dimensions 

When  implementing  upwind  weighting  in  two  dimen¬ 
sions,  several  additional  effects  occur.  Because  of  the 
multiplicative  generation  of  the  weighting  functions,  the 
y  direction  winding  affects  the  x  direction  equations 
substantially.  Thus,  when  considering  a  flow  with  a 


dominant  velocity  direction,  the  concept  of  cross  winding 
arises.  Velocity  in  this  cross  direction  significantly 
affects  the  response  of  the  main  flow. 

To  see  this,  partial  sums  of  terms  from  a  typical 
nodal  equation  will  be  presented.  The  elements  used 
will  be  bilinear  squares  with  unit  length  sides  (Fig.  4.2a). 
All  nodal  parameters  will  be  interpolated  linearly  and 
all  variables  will  use  the  same  interpolation  functions. 
Upwind  weighting  will  only  be  applied  to  the  convected 
terms . 

The  interpolation  functions  are 
N1  =  (1  -  x)  (1  -  y) 

N2  =  (x)  (1  -  y) 

N3  =  (x)  (y) 

N4  =  (1  -  x)  (y)  (4.2-6) 

The  weighting  functions  are 

w,  =  [0-x)  -  J«x0-idJ[0-yJ  - 
Wi-  [  x  -f  3*  i -*)][(  i- y)  -  Wy(i-y)] 

W5-f  x  +3*x(i -xj  y  +3fy(i- y)] 

Wy  s[0-x)  -5><x(i-x)][.  y  4  30y{/-y)J 


(4.2-7) 
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i  ■  "  "  ' 

l 

I  To  see  the  interaction  of  winding  in  the  two  direc- 

■  tions,  consider  the  nodal  equations  for  two  elements  with 
a  common  side  (Fig.  4.2b) .  The  eauation  for  v  at  node  3 

|  contains  convected  inertia  terms  that  are  quadratic  in  the 

1 

velocities  from  all  nodes.  For  purposes  of  comparison, 

•  I  2  2 

.  j  terms  in  v  ,v  v  ,v  from  the  v  and  v  equations 

i r  X  X  — .  X  j  X  a  X  a  X  a 

.  3  3  4  4  3  4 

will  be  listed. 

From  the  equation,  the  terms  are 

•  \ %  l d-W  f *  (i  7  + 

'^*7ps  t  (k ]  - 

.  +  \  Uk  1  *■ 

■  +(£-£»)  Lli  - 

|  '(hdl)  Ui  *-(k 

-({•'»$  [({-i?)  pt  +(i-i*)py]  + 

I  +  (s  f&*»«7  p,  *  (r *  **)  p»  ]* 

■  +&*»*)  P*l- 

'  -(i-iMUrb*) t,  f  li'i*)ps]- 

|  ft  +  (4  J /  + 

ft + 

^  *bf*)  fv  +  (&.*&  P,  ]  ~ 

-(&'»*) [(&-**) ft  Jj ■ 


I 

I 


(4 . 2-8) 


From  the  v  equation,  the  terms  are 
X4 

•  • .+  M*,  {&*&)[ +  (i  + 

+(&&)[ '  (i+h^Pil  - 
“i^ips  *  (h'hK)  fcl  ” 

+d’&r)fr  l}  + 

{(***$> (k+&*)p\  +(i*i*)ps1+- 

*(&*&*)[ Pv  +tj\+A>*)  Pzl  - 

-0>&)UU*)h  «A-b)9s1- 

-(^WC(4-tta  +  <*-Wfvl  + 

+(£*£*)  [^<W  p,  +  a  +h*)  ?z  7  + 
+(t*»9)[(i*l*)h  +6 U/Wpj  - 
M  'h^Pi  +(*•&*)  PrJ  " 

-  (i*iel)  [(h-to*)  &  +  (h  -If) pi  j( + 

*  V*  {  (&*•«)  [(i'Mfl  +(i*&*)p3l + 

4  (r  «»f)  I  (*  ♦£*)# *  fit  *iMf*  j  - 

**  "S^Ps  frJ  " 

[(/i’A")p4  +  -iWpyj/  f . . . 

(4.2-9) 


The  cross  winding  effect  can  be  seen  by  comparing 
the  terms  above  for  flow  in  the  x  direction  only.  For 
8  0,  the  terms  provide  unequal  contribution  to  the  v 

3 

and  v  equations.  A  one  row  strip  of  two-dimensional 
x4 

elements  cannot  solve  a  one-dimensional  problem  accurately 


unless  8=0. 
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To  further  illustrate  the  convected  inertia  terms 
in  the  typical  equation,  selected  terms  from  the  v 

y3 

equation  will  be  presented.  The  particular  terms  involv¬ 
ing  the  y  velocities  are  given,  followed  by  simplifica¬ 
tions  to  clarify  their  meaning. 

Let  all  nodal  values  of  the  density  take  on  the 
uniform  value  c q-  First  consider  a  =  6  =  0,  or  no  upwind 

weighting.  With  these  simplifications,  the  v  terms 

y3 

reduce  to 


ytt  {(2v  +v  +3v  )  (v  -v  )  +(v  +2v  -12v  -6v  + 

144  y-^  y  2  Y3  y2  yl  yl  y2  y3  y4 

+  v  +  2v  )  (v  -v  )  +  (3v  +2v  +v  )  (v  -v  )} 

y5  y6  y4  y  3  y3  y5  y6  y6  y3 

(4.2-10) 

Notice  that  the  appearance  is  that  of  weighted  sums  of 
the  three  available  velocity  gradients. 

Letting  S  =  1  and  a  =  0  corresponds  to  a  dominant  y 
velocity  across  the  element  row.  The  terms  then  reduce  to 


ytt  {(v  -v  +v  -v  ) (v  -v  )+(v  -v  +6v  -6v  + 

144  y1  y  2  y3  y4  y2  y1  yx  y2  y3  y4 


+v ,  -v  )  (v  -v  }  +  {V  -V  +V  -V  )  (v  -V  )} 

y5  y6  y4  ^3  y 3  y4  y5  y6  y6  y5 

(4.2-11) 


Although  the  form  of  weighted  sums  of  velocity  gradients 

can  still  be  seen,  the  weights  have  changed.  Terms  on  the 

upwind  row,  v  ,v  , v  ,  are  slightly  increased,  while 
yl  y3  y5 
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terms  on  the  downwind  row, 


are  substantially 


reduced . 

Finally,  consider  a  =  B  =  1  or  dominant  velocities 
in  the  positive  direction  of  each  axis.  The  terms  reduce 
to 


-s -j-%  { ( 14v  -  14v  )  { v 

144  y.  y,  y. 


VV  >  + 
yl 


+  (  -22v  +22v  -30v  -2v  +2v  )  (v  -v  )  + 

yl  y  2  y4  y  5  y6  y4  y  3 


+  28v  v  +  ( -  4v  +4v  )  (v  -v  P  (4.2-12) 

y 3  y 4  y5  y6  y6  y5 


Here  the  shift  to  upwind  is  less  obvious. 


Upwind  Weighting  on  Terms  other  than  Convected  Inertias 
In  a  consistent  application  of  the  present  varia¬ 
tional  method,  the  same  weighting  functions  are  used  in 
all  terms  of  an  equation.  To  assure  stability,  however, 
it  appears  necessary  only  to  apply  upwind  weighting  to  the 
convected  terms.  In  practice,  nodal  equations  of  elements 
without  an  upwind  neighbor  are  ill-conditioned  by  the 
application  of  such  weighting  to  the  entire  equation. 

This  difficulty  may  be  avoided  by  prescribing  nodal 
values  for  upwind  boundary  nodes.  Such  artificial  bound¬ 
ary  conditions  are  not  desired  since  it  may  be  difficult 
to  find  appropriate  values  to  prescribe  prior  to  solu¬ 
tion  of  the  problem. 
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To  illustrate  this,  elements  of  the  density  rate 
matrix  [RM]g  will  be  presented,  where 

[RM]  =  /  [lV[L]  dV  (4.2-13) 

F 

e 

Using  the  bilinear  interpolation  functions  for  [L]  and 

★ 

the  upwinded  weighting  functions  for  [L  ] ,  the  elements 
of  the  first  row  are 

RM (1,1)  =  (i  -  ~  a)  (i  -  i  8) 

RM  (1,2)  =  (|  -  |  a)  (j  -  |  6) 

RMd,3)  =  (|  -  i  a)  (£  -  j  8) 

RM  (1,4)  =  (i-  -  i  a)  (|  -  |  8)  (4.2-14) 

Consider  the  case  of  full  a  winding,  and  no  8  winding. 

Two  of  the  terms  become  negative,  while  following  a 
lumping  procedure,  the  sum  of  these  terms  is  zero.  When 
assembled  into  a  global  matrix,  there  will  be  no  nonzero 
terms  added  if  there  are  no  upwind  neighbors. 

The  momentum  rate  or  mass  matrix  behaves  similarly. 
This  matrix  is 

[VM]  =  f  [N*]T[N]  [L]  {p}  dV  (4.2-15) 

Fe 

Using  the  bilinear  functions  for  [N]  and  [ L ] ,  and  the 

* 

winded  functions  for  [N  ] ,  the  lumped  mass  is 


V 
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VM  (1,1)  - - —  p  (4.2-16) 

360 

where  p  is  an  assumed  uniform  value  of  the  density, 
and  a  =  1,  8=0.  For  comparison,  the  result  with  no 
winding  and  uniform  density  is  — p  for  the  lumped  mass. 


4.3  Discretization  of  the  Elastic  Solid  Functional 
Membrane  Element 

The  elemental  contributions  arising  from  the 
partial  differentiation  of  the  leading  term  in  the  solid 
functional  with  respect  to  the  nodal  values  of  the  adjoint 
displacement  will  be  called  the  element  restoring  forces. 
This  membrane  element  will  have  structural  displacements 
that  are  linear  combinations  of  the  two  sets  of  nodal 


The  adjoint  variables  will  be  interpolated  similarly. 

Evaluating  the  strains  and  assuming  a  linear  elastic 
material,  the  stresses  are  found  to  be 


’■v 
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s  =  Ee 

XX  XX 


s  =  s  =  2Ge 
xy  yx  xy 

s  =0 

yy 


(4.3-2) 


The  leading  term  in  the  element  functional  is  found, 
after  insertion  of  the  precedina,  to  be 

4 


=  [  *4,*  4  ^ f  +  **,yJ  ^  ^ 


i 


'St 

t>  i 


C 


«„>+ 


4  (  uy,  +  S,x  *  *Ky)J  h 


The  element  restoring  forces  are 


(4.3-3) 


a< 


K 


L 

=  Q*, =  /  i  h ^ ^ V +  ^  i:  uv,x + 

+•  N,*  ^  +  2  +  I  WY.a  ^x,*  ^  (4 . 3-4a) 

~  Qy  ~  ■/  fw.,»  ^  +  ^l,A  ^  ^Y,X  + 

+  i  <,  U*,  +  i  )]  <U 


( 4 . 3-4b) 
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4N,‘.hefu'.  ti 


(4 . 3-4c) 


+  Nj*  ^  ^  ^  u*,x  Uy(  x 

4**)^  ( 4 . 3-4d) 


These  results  are  identical  to  those  obtained  by 
retaining  all  terms  in  the  virtual  work  equations  for 
large  deformations.  To  reach  the  simplified  form  usually 
used  in  practice,  the  assumptions  required  are 

1  .  1  +  u  =1 

X  ,X 

2-  \  ux,x  c  '  UX,X 

3-  sxy  =  0  but  uy/X  +  0 

These  assumptions  are  consistent  with  the  practice  of 
neglecting  higher  powers  of  small  terms. 
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4.4  Fluid  Structure  Interface 

Typical  Nodal  Equation  on  the  Interface 
The  condition  on  the  interface  S^.  of  the  fluid 
region  Vp  and  the  structural  region  Vg  is  continuity  of 
velocity  and  temperature.  For  the  specific  case  of 
inviscid  flow,  this  is  reduced  to  continuity  of  the 
normal  velocity.  In  practice,  this  means  the  equation 
for  the  normal  velocity  of  a  node  on  the  interface  will 
contain  contributions  from  elements  in  both  regions. 

The  example  presented  here  is  the  momentum  equation 
for  a  node  common  to  two  unit  square  fluid  elements  and  two 
membrane  elements  (Fig.  4.3).  This  equation  will  contain  con- 

vected  inertia  and  pressure  terms  from  the  fluid  element 
and  restoring  forces  from  the  membrane.  In  addition,  the 
nodal  mass  will  be  the  sum  of  lumped  masses  from  the 
fluid  and  membrane  elements. 

The  equation  for  the  y  velocity  at  node  3  will  be 
presented.  An  external  pressure  of  magnitude 
p  will  be  applied  to  the  membrane.  Each  membrane  element 
will  have  a  thickness  h  and  uniform  density  pg.  In  this  exam¬ 
ple,  contributions  from  the  x  velocities  to  the  convected 
inertias  will  be  neglected.  The  fluid  will  have  a  uniform 
density  and  pressure  p^.  These  conditions  are 
appropriate  for  initial  motion  of  the  membrane  into 
stationary  fluid. 

Some  contributions  to  the  equation  for  v 


are : 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


structural  mass  = 
fluid  lumped  mass 


pressure  gradient 


2  2  PS^1  psh 

1  ^  1  ,  2  .1 
=  T8  pl+  36  °2+  9  p3  +  9  p4  + 

1  1  _  1 
18  p5  36  P6  2  PF 

JL  _  _1_  1  1 

~  12  P1  ~  12  P2  ~  3  P3  ”  3  P4  ~ 

12  P5  “  12  P6  =  '  P0 


external  pressure 


Let  a  =  0 ,  £  =  1,  and  Q  be  the  spring  restoring  force. 

y3 

The  equation  for  v^  is 

(psMp,JV  “‘Vi0- +  P  ~ 
“SfctVVVV'VV  +  (W4V 


(4.4-1) 


If  p  has  the  value  of  the  rest  pressure  of  the  fluid, 
there  is  no  pressure  gradient  across  the  membrane.  If  the 
fluid  velocities  away  from  the  membrane  can  be  neglected, 
the  equation  simplifies  to 


(fth  9f\  •  -vl  M  J 


(4.4-2) 
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I 

I 

I 

I 

I 


I 

I 

I 


f 


Continuity  on  the  Interface 

The  only  condition  on  admissible  trial  functions 

on  the  interface  surface  S^.  is  that  the  velocity  and 

temperature  fields  be  continuous.  Continuity  of  trial 

functions  is  a  common  requirement  in  finite  element 

procedures  and  is  generally  met  by  requiring  elements 

to  be  compatible.  The  following  will  investigate  the 

compatibility  of  a  structural  element  and  a  fluid  element. 

To  facilitate  this,  consider  the  mesh  segment 

shown  in  Fig.  4.4.  The  elements  2,3,4  and  8  are  fluid 

elements,  elements  1,5,6  and  7  are  structural  elements, 

and  nodes  1,5  and  9  are  the  interface  surface  nodes. 

Let  uq  represent  the  structural  displacements,  va  be 

the  fluid  velocities,  be  the  interpolation  function 

for  node  i  in  element  e,  and  Ne  denote  that  function 

1S 

restricted  to  the  interface.  The  triangular  elements 
will  allow  linear  variation  in  the  variable  fields  and 


this  will  be  assumed  for  both  regions.  Let  u_  be  the 


function  u 


restricted  to  element  e  on  surface  S, 


Then  in  element  8,  a  fluid  element,  the  velocity  is 


v  =  Z  N8  v 

a  .  i 


i  =  5,6,9 


=  N?  v  +  N?v  +  N? v 
5  a,  6a  9  ag 
3  6  ~ 


(4.4-3) 


0 

If  N  =0  as  is  common,  the  velocity  on  S  is 
6S  1 


v8  =  N8  v  +  N®  v 
aS  5S  a5  9S  ap 


(4.4-4) 
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In  a  similar  manner,  the  temperature  is 


Tg  =  M®  T5  +  M®  Tg 
S  S 


(4.4-5) 


In  element  7,  a  structural  element,  the  displacement  is 


u7  =  N?u  +  nZu  +  nZu 
a  5  a5  8  cig  9  ag 


(4.4-6) 


Then,  on  Sj,  the  displacement  becomes 


7  7  7 

u  =  Nc  u  +  N'u 

aS  5S  a5  9  a9 


(4.4-7) 


N8 

s 


and  the  temperature  is 


7  7  7 

TS  =  M5  T5  +  M9  T9 
S  S 


(4.4-8) 


Since  Ne  =  0  at  node  j  f  i  and  Ne  =  1  for  node  i,  at 
i  J  l 

node  5,  the  interface  values  of  the  field  variables  are 


v  =  v 
aS  a5 

T8  =  T 
lS  5 


(4.4-9) 


7 

u  =  u 
aS  a5 

T7  =  T_ 


(4.4-10) 


Since  the  interpolation  functions  are  not  time  dependent, 


the  structural  velocity  is 
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elements . 
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4.5  Time  Integration 

The  actual  form  of  the  matrix  equations  is  an 
integral  over  time,  e.g., 

— -  V"  =  f+2  f  [VMJ  (v)  +  [WT  iv)  -{P)1  dt  =  0 

3  {v  tl  L  J  (4.5-1) 

This  integral  equation  is  satisfied  in  practice  by 
solving  the  ordinary  differential  equation  obtained  by 
setting  the  integrand  to  zero  for  all  time. 

For  reasons  of  simplicity  in  implementation,  a 
Runge-Kutta  integration  package  was  used.  This  routine 
calculates  end  of  time  step  values  based  upon  a  weighted 
average  of  the  rates  obtained  from  four  evaluations 
during  the  step. 

Although  this  is  a  highly  nonlinear  equation  set,  the 
time  step  criterion  for  linear  problems  has  proven 
adequate.  Specifically,  the  time  step  is  chosen  such 
that  information  propagating  at  the  local  speed  of  sound 
travels  no  farther  than  one  element  width  in  one  time 
step.  In  practice,  one  sound  speed  is  chosen,  usually 
the  rest  speed  of  sound.  The  maximum  time  step  is 
calculated  based  upon  the  smallest  element  in  the  field. 
Since  local  sound  speeds  may  rise  considerably  above  the 
rest  value,  the  calculated  time  step  is  halved. 

Parametric  studies  of  numerical  solutions  have  shown 
this  to  be  adequate. 


55 


5.  NUMERICAL  EXAMPLES 

Several  numerical  problems  are  presented.  Their 
purpose  is  to  exercise  the  major  capabilities  of  the 
variational  formulation  and  its  computer  implementation. 
Since  a  combined  problem  with  compressible  isentropic 
fluid  and  a  large  deformation  structure  is  quite  complex, 
simpler  examples  are  given  first.  Understanding  of  these 
simple  examples  was  essential  to  adequately  solve  the 
combined  problem. 

To  verify  the  fluid  element  used  in  the  present 
method,  problems  of  a  one-dimensional  wave  tube  have 
been  solved.  This  problem  can  demonstrate  extreme  fluid 
motions  for  which  a  limited  analytical  solution  is 
available . 

The  wave  tube  is  enlarged  to  two  dimensions  by 
examining  a  radially  expanding  cylinder  in  an  infinite 
fluid.  Valuable  insight  into  the  functioning  of  the 
upwind  weighting  is  obtained  here. 

The  final  example  is  transient  motion  of  a 
membrane  with  fluid  on  one  side.  In  the  first  of  two 
cases  considered,  initial  motion  of  th^  membrane  dis¬ 
turbs  the  fluid.  In  the  second,  a  pressure  wave 
traverses  the  fluid  and  excites  the  membrane. 


56 


5.1  Piston  Wave  Tube 

The  one-dimensional  gas  dynamics  of  a  rigid  piston 
moving  in  a  tube  filled  with  a  perfect  gas  was  investigated. 

The  pressure  and  velocity  of  the  inviscid  gas  were  calcu¬ 
lated  using  a  prescribed  motion  of  the  piston.  The  move¬ 
ment  of  the  piston  is  modeled  by  changing  the  coordinates 
and  velocities  of  the  nodes  assigned  to  the  piston 
surface.  When  the  element  in  front  of  the  piston  is 
reduced  to  a  certain  size,  it  is  eliminated  from  the 
problem  solution. 

The  rest  speed  of  sound  for  the  following  examples 
is  1000  fps.  In  the  first  example,  the  piston  is  started 
impulsively  with  a  speed  of  500  fps  and  keeps  this 
speed  for  all  time.  I  i  the  second  example,  the  piston 
velocity  increases  linearly  from  rest  to  500  fps  at  .002 
sec,  then  decreases  linearly  to  rest  at  .004  sec.  The 
parameter  a  is  the  upwind  weight  in  the  direction  of 
motion.  With  one  exception  which  will  be  noted,  symmetric 
weighting  was  used  in  the  transverse  direction.  The  time 
step  was  half  the  linear  stability  limit  based  upon  the 
rest  speed  of  sound  and  the  initial  configuration  mesh 
size. 

Figure  5.1  presents  the  fluid  velocity  and  pressure 
for  the  constant  speed  piston.  The  calculated  pressure 
ratio  across  the  shock  agrees  to  less  than  1%  with  that 

V 
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given  by  tt\£  Rankine-Hugoniot  relations  (Liepman,  Roshko, 
1957).  The  smearing  of  the  shock  discontinuity  is 
typical  of  the  finite  element  solutions.  There  appears 
to  be  very  little  ringing  or  oscillation  in  the  solution 
behind  the  shock. 

The  results  for  the  second  piston  example  are 
presented  in  Fig.  5.2.  The  example  contains  a  compres¬ 
sion  wave  which  changes  to  a  shock  as  well  as  an  expan¬ 
sion  wave.  Also  presented  are  some  early  time  solutions 
by  the  method  of  characteristics  for  comparison. 

Agreement  between  the  two  is  quite  good  in  the  expansion 
wave  but  due  to  smearing  of  the  finite  elements, 
agrees  only  fairly  in  the  high  gradient  region  of  the 
shock  front.  Again  there  is  little  ringing  near  the 
peak,  and  the  fluid  behind  the  wave  is  reasonably  still. 

The  deterioration  of  the  solution  as  a  is  decreased 
to  zero  is  presented  in  the  next  figure  (Fig.  5.3). 
Deviations  from  the  full  upwind  solution  (a  =  1.0)  are 
noticeable  in  the  case  of  a  =  .5.  The  solution  for  a  *  0 
is  barely  recognizable.  These  results  are  consistent  with 
those  of  Zienkiewicz  (1977) . 

It  is  a  characteristic  of  the  present  upwind 
weighting  scheme  that  winding  in  one  direction  affects 
the  solution  in  the  other  direction.  An  example  of  this 
is  presented  in  Fig.  5.4.  The  calculated  velocities  for 
6  *  ±  1  differ  significantly  from  the  6=0  solution 
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of  Fig.  5>  2  . 

5.2.  Cylindrical  Wave 

To  demonstrate  the  two-dimensional  capability  of 
the  present  method,  the  piston  problem  was  expanded  to  a 
cylindrical  problem  (Fig.  5.5).  The  goal  is  to  calculate 
velocities  and  pressure  in  an  initially  still  fluid 
surrounding  an  expanding  cylinder.  The  cylinder  radius 
and  velocity  are  prescribed  functions  of  time. 

The  init-  1  radius  of  the  cylinder  is  1.0  ft.  Its 
velocity  increases  linearly  to  500  fps  at  .002  sec  then 
decreases  linearly  to  zero  at  .004  sec.  The  rest  speed 
of  sound  is  1000  fps.  The  physical  problem  is  in  fact 
one-dimensional  in  the  radius.  A  solution  procedure 
using  the  method  of  characteristics  is  given  by  Rudinger 
(1969)  but  was  not  used  because  of  its  complexity. 

Figure  5.5  also  presents  the  calculated  radial  veloc¬ 
ity  and  pressure.  Since  there  should  be  no  tangential 
velocity,  the  appropriate  value  of  the  tangential  or  cross 
wind  weighting  is  zero.  Table  5.1  shows  the  effect  on 
the  main  flow  caused  by  the  inclusion  of  cross  wind 
weighting.  The  local  value  of  the  winding  parameter  is 

6  =  1 . 0  x  sign  (v  ) 


The  variation  in  the  radial  velocity  and  pressure  are 
evidently  associated  with  the  fluctuating  sign  of  the 
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very  small  tangential  velocity. 

5.3  Transient  Motion  of  a  Membrane 

The  most  complex  problem  attempted  was  transient 
motion  of  a  membrane  wetted  on  one  side  by  an  inviscid 
compressible  fluid.  Such  a  problem  exercises  all  the 
features  of  the  computer  program  and  illustrates  the 
problems  that  can  be  attacked  by  the  present  method.  The 
membrane  is  capable  of  accurately  reproducing  large 
motions  and  the  fluid  permits  accurate  and  stable  calcu¬ 
lation  of  motion  in  the  presence  of  extreme  gradients. 

The  problem  at  hand  will  illustrate  the  proper  coupling 
of  the  two  regions  at  their  interface. 

The  membrane  is  stretched  at  y  =  0  between 
-2.5  £  x  2 . 5  with  rigid  extensions  to -5.0  and  7.5  in 
the  x  direction  (Fig.  5.0).  The  fluid  volume  is  truncated 
at  y  =  5.0  and  extends  from  x  =  -  5.0  to  x  =  7.5,  The 
membrane  starts  with  no  deflection  and  an  initial  velocity 
such  that  a  linear  membrane  in  a  vacuum  would  have  a 
maximum  deflection  of  8%  of  its  length.  The  fluid  will 
start  either  at  rest  or  have  a  uniform  velocity  in  the  x 
direction  of  magnitude  v^  =  1000  fps  or  M^  =.31.  The 
rest  sound  speed,  a^,  is  3200  fps,  and  the  rest  pressure, 
Pq,  is  2000  psf.  The  membrane  has  a  small  amplitude  in 
vacuum  natural  frequency  of  800  cps.  A  uniform  pressure 
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of  magnitude  equal  to  the  fluid  rest  pressure  is  applied 
to  the  outside  of  the  membrane. 

Table  5.2  presents  variables  at  the  membrane  at  a 
fixed  time.  The  membrane  has  moved  into  the  fluid  to  a 
point  just  short  of  its  maximum.  In  all  cases,  the 
membrane  is  never  affected  by  the  presence  or  absence  of 
the  free  stream  velocity.  The  fluid  velocity  and  pressure 
show  a  marked  asymmetry  when  the  cross  wind  is  present. 
Apparently  the  asymmetric  pressure  gradient  and  x  velocity 
are  too  small  to  affect  the  membrane  motion  in  so  short 
a  time. 

The  physical  properties  of  the  present  example  of 
fluid  and  structure  have  been  chosen  such  that  the  membrane 
is  not  appreciably  affected  by  the  presence  of  the  fluid. 
For  example,  the  period  of  oscillation  is  shortened  by 
only  a  few  percent,  and  the  amplitude  is  not  affected. 

Since  the  fluid  has  been  truncated  and  a  prescribed 
normal  velocity  imposed,  reflections  from  these  boundaries 
will  lead  to  a  solution  that  looks  more  like  vibrations 
in  a  cavity  at  the  later  times.  Some  effects  attributable 
to  reflections  can  be  seen  in  the  results  that  follow. 

Figure  5.7  presents  pressure  wave  forms  along  the 
cavity  center  line.  The  times  chosen  correspond  to 
quarter  periods  of  the  membrane  motion.  The  disturbances 
appear  to  be  large  and  the  distortion  in  the  wave  shape 
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can  be  attributed  to  the  changing  local  speed  of  sound. 

The  amplitudes  of  the  third  cycle  (compare  t  =  2.76  x  10  3 

-3  -3 

to  t  =  1.50  x  10  and  t  =  0.30  x  10  )  are  noticeably 

larger  than  the  first.  This  is  probably  due  to  the 
reflections  off  the  boundary  at  y  =  5.0. 

Pressure  wave  forms  on  two  planes  parallel  to  the 
membrane  show  similar  effects  of  reflections  at  later 
times  (Fig.  5.8).  In  addition,  results  of  a  calculation 
with  an  initial  velocity  in  the  x  direction  of  1000  fps  or 

M  =  .31  have  been  included.  At  all  times  presented, 

00 

the  waves  have  drifted  noticeably  downstream. 

As  a  second  case  using  the  same  model  geometry, 
response  of  the  membrane  to  a  traveling  pressure  wave  was 
considered.  Both  fluid  and  membrane  are  initially  at  rest 
and  a  step  pressure  increase  from  2000  psf  to  4000  psf  is 
applied  to  the  fluid  boundary  at  y=5.0.  The  wave  propa¬ 
gates  across  the  cavity  and  reflects  off  the  membrane  and 
its  rigid  extensions. 

Three  typical  pressure  wave  forms  on  the  cavity 
center  line  are  presented  in  Fig.  5.9.  The  times  chosen 
depict  the  incident  wave,  the  time  of  maximum  pressure  and 
the  wave  as  it  returns  to  the  y=5.0  boundary.  The  motion 
of  the  center  of  the  membrane  is  shown  in  Fig.  5.10.  In 
Table  5.3  this  same  information  is  presented  for  compar¬ 
ison  to  a  solution  obtained  by  modal  expansion.  For  the 
comparison  solution,  the  motion  of  the  membrane  was  calcu¬ 
lated  using  the  pressure  at  y=0  as  the  forcing  input. 
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LIST  OF  SYMBOLS 

A  free  energy 

specific  heat  at  constant  volume 
E  internal  energy  of  the  fluid  or  Young's  modulus 
G  shear  modulus 
I  functional 

L  density  interpolation  function 

M  Mach  number  or  temperature  interpolation  function 
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U  internal  energy  of  the  solid 
W  weighting  function 

a  local  sound  speed  or  acceleration  vector 
d  spatial  rate  of  strain  tensor 
e  material  strain  tensor 
f  body  force  vector 
g  metric  tensor 

h  membrane  thickness  or  mesh  spacing 
k  thermal  conductivity 
n  surface  normal 
p  pressure 

q  velocity  magnitude  or  heat  flux  vector 
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r  rate  of  heat  generation  or  radial  coordinate 
s  material  stress  tensor 
u  displacement 
v  velocity 

a  £  or  x  direction  winding  parameter,  ~  1  £  a  £  1 
£  n  or  y  direction  winding  parameter,  -  1  _<  £  <_  1 
Y  gas  constant 
6  Kronecker ' s  delta 

a  viscosity  coefficient  or  Lame  constant 
u  viscosity  coefficient  or  Lame  constant 
£  coordinate  of  quadrilateral  element 
n  coordinate  of  quadrilateral  elemer  t 
o  density 

$  velocity  potential  or  direction  in  Frechet  differential 
ip  direction  in  Frechet  differential 
e  cylindrical  coordinate  angle 

Superscripts  and  Subscripts 
S  solid 
F  fluid 

a,b,c,  .  . .  indices  of  material  coordinates 
e  element  quantity 
r  radial  component 
t  tangential  component 

indices  of  rpatial  coordinates 

oc  far  field  value 

* 

f  adjoint  quantity  f 


f  prescribed  quantity  f 


*4H  4J 

P  P 
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Operators 

Vf  gradient  of  f 

f>x  partial  differentiation  of  f  with  respect  to  x 
f|x  covariant  differentiation  of  f  with  respect  to  x 

material  derivative  of  f 

f  time  derivative  of  f 
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2251 
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1.2 
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22.5 
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-0.2 

2700 
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0.3 
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-1.9 
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2254 

308.4 

0.5 

2483 

67.5 

374.0 

0.2 

2701 

307.5 

-0.3 

2480 

78.5 

234.9 

-3.8 

2251 

307.5 

-1.2 

2479 

90.0 

384.9 

0.0 

2733 

307.8 

0.0 
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Note:  p0 

=  1700  psf,  ag 

=  1000 

fps  and 

Ar  =  . 

125  ft 

Table  5.1.  Cylindrical  wave  solution  at  t  ■  .003  sec 


r  *  3.0  ft,  with  and  without  cross  wind  weighting 
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-2.0 
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-1.5 

.2305 
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.3179 

886 

-178 

2540 

2218 

-0.5 

.3744 
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0. 
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0. 
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0. 

2000 

2000 

Note:  Pg  ■ 

2000  psf 
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• 

Table  5.2.  Fluid  velocity  and  pressure  on  membrane 
at  t  *  .30  x  10‘3. 
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Time 

2  (sec) 

Displacement 

Finite  Element 

xlO3  (ft) 

Displacement 
Modal  Expansion 
xlO3  (ft) 

132 

-.258 

-.263 

144 

-.573 

-.601 

156 

-.993 

-1.08 

168 

-1.40 

-1.58 

180 

-1.62 

-1.88 

.192 

-1.56 

-1.81 

,204 

-1.17 

-1.34 

,216 

-.599 

-.633 

.228 

-.058 

.030 

.240 

.237 

.385 

.252 

.184 

.304 

.264 

-.183 

-.159 

.276 

-.728 

-.814 

.288 

-1.24 

-1.42 

.300 

-1.54 

-1.76 

Table  5.3.  Membrane  response  to  incident  pressure  wave 
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